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Abstract 



Context. Mass loss plays a dominant role in the evolution of low mass stars while they are on the Asymptotic Giant Branch (AGB). The gas 
and dust ejected during this phase are a major source in the mass budget of the interstellar medium. Recent studies have pointed towards the 
J—, importance of variations in the mass-loss history of such objects. 
Or Aims. By modelling the full line profile of low excitation CO lines emitted in the circumstellar envelope, we can study the mass-loss history of 
q . AGB stars. 

?H Methods. We have developed a non-LTE radiative transfer code, which calculates the velocity structure and gas kinetic temperature of the 
jy-j envelope in a self-consistent way. The resulting structure of the envelope provides the input for the molecular line radiative calculations 
& • which are evaluated in the comoving frame. The code allows for the implementation of modulations in the mass-loss rate. This code has been 
^ ' benchmarked against other radiative transfer codes and is shown to perform well and efficiently. 
. , Results. We illustrate the effects of varying mass-loss rates in case of a superwind phase. The model is applied to the well-studied case of VY 
■ CMa. We show that both the observed integrated line strengths as the spectral structure present in the observed line profiles, unambiguously 
5— l demonstrate that this source underwent a phase of high mass loss (~ 3.2 x 10~ 4 M Q yr -1 ) some 1000 yr ago. This phase took place for some 
100 yr, and was preceded by a low mass-loss phase (~ 1 x 10~ 6 M Q yr ) taking some 800 yr. The current mass-loss rate is estimated to be in 
the order of 8 x 10~ 5 M Q yr~' . 

Conclusions. In this paper, we demonstrate that both the relative strength of the CO rotational line profiles and the (non)-occurrence of spectral 
structure in the profile offer strong diagnostics to pinpoint the mass-loss history. 

Key words. Line: profiles, Radiative transfer, Stars: AGB and post-AGB, (Stars): circumstellar matter, Stars: mass loss, Stars: individual: VY 
CMa 



1. Introduction 

When low and intermediate mass stars approach the end of 
their lives and ascend the Asymptotic Giant Branch (AGB), 
mass loss dominates the subsequent evolution ultimately lead- 
ing to the removal of the entire envelope. In this process, 
the star develops a dust-driven wind with a velocity of 10 - 
15kms _1 (Habing & Olofsson 2003, and references therein). 
AGB stars lose ~ 35 - 85 per cent of their mass through the 
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stellar wind before ending up as a white dwarf (Marshall et al. 
2004). Likewise, more massive stars (Mzams ^ 8 M ) may 
pass through a red supergiant phase and lose mass in a simi- 
lar manner. 

The mass-loss properties change with time as these stars 
ascend the AGB and become more luminous. It is often pos- 
tulated and in rare cases observed (Justtanont et al. 1996; van 
Loon et al. 2003; Marshall et al. 2004, e.g.) that the AGB evo- 
lution ends in a very high mass-loss phase, the so-called su- 
perwind phase (Iben & Renzini 1983). The nature and onset of 
such a superwind are far from understood. 

This superwind phase is not the only type of variable mass 
loss. Recent observations of (post) AGB objects and Planetary 
Nebulae (PNe) show that winds from late type stars are far from 
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being smooth. The shell structures found around e.g. CRL 2688 
(the Egg Nebula, Ney et al. 1975; Sahai et al. 1998), NGC 6543 
(the Cat's Eye Nebula, Harrington & Borkowski 1994) and the 
AGB star IRC +10 216 (Mauron & Huggins 1999, 2000; Fong 
et al. 2003), indicate that the outflow has quasi-periodic os- 
cillations, with density contrasts corresponding to mass-loss 
variations of up to a factor -100-1000 (Schoier et al. 2005). 
The time scale for these oscillations is typically a few hundred 
years, i.e. too long to be a result of stellar pulsation, which 
has a period of a few hundred days, and too short to be due 
to nuclear thermal pulses, which occur once in ten thousand 
to hundred thousand years. Presently, the physical mechanism 
responsible for these variations, their rate of occurrence and 
their importance in terms of the amount of ejected matter in- 
volved are unknown. Possibly they are linked to the hydrody- 
namical properties of a dust-driven wind. For example, Simis 
et al. (2001) find that a feedback mechanism between the wind- 
acceleration and dust growth at the base of the wind may cause 
modulations in the mass-loss rate. 



Already in the late seventies, low excitation rotational 
carbon-monoxide (CO) emission lines emerging from the cir- 
cumstellar envelopes (CSEs) around AGBs were observed and 
analysed to estimate the mass-loss rates during the AGB phase 
(e.g. Zuckermann et al. 1977; Morris 1980; Knapp et al. 1980, 
1982; Knapp & Morris 1985). CO has proven to be instrumen- 
tal for this since (almost) all the elemental carbon is locked up 
in CO throughout most of the envelope. While previous studies 
have assumed a constant mass-loss rate in the analysis of rota- 
tional CO lines, the goal of our study is to exploit the rotational 
CO line profiles of different excitation levels as a complemen- 
tary diagnostic to probe variations in mass-loss rate. With the 
arrival of HIFI and PACS the powerful spectrometers on the 
Herschel Space Observatory — ESA's far infrared mission to 
be launched in 2008 — and ALMA — the submillimeter in- 
terferometer which is being built by ESO and NRAO in Chile 
— new high-quality data of CO emission lines will put new 
light on this ongoing research of understanding the geometry, 
physical conditions and mass-loss history of AGB stars. 



In this paper, we describe our non-LTE (non- Local 
Thermodynamic Equilibrium) radiative transfer code devel- 
oped to predict the CO molecular line strengths. In Sect. 2, the 
equations governing the temperature and velocity structure in 
the CSE are presented. Since the CO rotational line profiles 
are very sensitive to the temperature structure in the CSE, the 
kinetic temperature is calculated self-consistently taking the 
main heating and cooling processes into account (Sect. 2.2). 
The code can treat any density profile, allowing us to deduce 
the amplitude of the mass-loss variations (Sect. 2.4). The non- 
LTE radiative transfer code is described in Sect. 3, and the ef- 
fect of variations in mass loss on the predicted CO profiles is 
studied. In Sect. 4, we apply our theoretical code to model the 
observed CO rotational emission lines in the supergiant VY 
CMa. We summarise the results in Sect. 5. 



2. Theoretical model for the temperature and 
velocity structure in the CSE 

A detailed understanding of the wind around red giants requires 
the solution of the equations of hydrodynamics, specifically in 
the complex region starting from the stellar surface into the re- 
gion of dust formation. The equations should include the con- 
servation laws of mass, momentum, and energy, in combina- 
tion with a chemical model describing the gas phase chemistry 
in the CSE. In the absence of such a theory, we follow previ- 
ous work, and focus on those equations describing the velocity 
and temperature structure in the CSE. The main assumption 
is a spherically symmetric mass loss. The code builds on the 
work of Justtanont et al. (1994, hereafter referred to as J94). 
Main changes are the inclusion of new and improved heating 
and cooling rates (Sect. 2.2), and the implementation of a vari- 
able mass loss (Sect. 2.4). Helium, in addition to atomic and 
molecular hydrogen, is now also included as collision partner. 

2.1. The gas and dust velocity 

In this study, we will focus on red giants with outflow ve- 
locities in excess of 5kms _1 , and mass-loss rates larger than 
3x 10" 7 M yr _1 . These winds are believed to be driven by radi- 
ation pressure on the dust which condenses in the outer parts of 
the atmosphere (Winters et al. 2003). In this situation, the equa- 
tions expressing the conservation of mass and momentum are 
respectively given by (e.g. Goldreich & Scoville 1976; Tielens 
1983, hereafter referred to as GS76 and T83 respectively) 

=M(r) = 4nr 2 p(r)v(r), (1) 

and 

,Jv{r) GM* 
v(r)— — = (F(r) - 1) — — , (2) 
dr r 1 

where M(r) refers to the gas mass-loss rate at a radial distance 
r from the star, p(r) is the gas density, v(r) the gas velocity, 
M* the stellar mass, and F(r) the ratio of the radiation pressure 
force on the dust to the gravitational force. 

The total dust mass-loss rate, M d , can be written as (J94) 

M d (r) = j ' M d {a,r)da 

= J A(r)a- X5 n H (r)j7Ta 3 pAxr z [v(r) + v A nf t (a,r)]da, (3) 

where Md(a, r) represents the dust mass-loss rate of a grain of 
size a, n H the hydrogen number density, p s the specific den- 
sity of dust and v^ia, r) the drift velocity of a grain of size 
a. For the dust composition, we use silicate grains with p s = 
3.3gcrrT 3 (Draine & Lee 1984). We adopt a grain-size distri- 
bution n d (a, r) da = A(r) cT 3 - 5 n H (r) da, with a minimum grain- 
size, fl m i„, of 0.005 pm and a maximum grain-size, a max , of 
0.25 pm. The slope -3.5 is typical for interstellar dust (Mathis 
et al. 1977). It is unclear whether this slope is representative 
for circumstellar dust shells l . Studies on dust condensation 

1 The choice of slope has very little direct influence on the results. 
The calculated gas-to-dust ratio changes to some extent since the ab- 
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and growth suggest that the circumstellar dust grain size distri- 
bution could be much steeper (Dominik et al. 1989). However, 
shattering by grain-grain collisions could lead to a much shal- 
lower power law (Biermann & Harwit 1980). The quantity A(r) 
represents an abundance scale factor giving the number of dust 
particles in units of particles per H atom. For the ISM, A is a 
constant, estimated to be 10~ 2510 cm 3 5 /H for silicate grains and 
2Q-25.13 cm 3-5^j_[ f or car b on grains (Draine & Lee 1984). In the 
case of CSEs around giants, the value of A(r) may differ from 
the ISM value, since it incorporates the (unknown) dust-to-gas 
mass ratio (if/) of the CSE. Moreover A(r) varies as a function 
of the radius r approximately as v(r)/(v(r) + Vdrift(O)' In case 
of a constant mass-loss rate, A decreases rapidly at the base 
of the envelope w.r.t. the value when the dust first condenses 
and the dust velocity is still zero, as the dust is rapidly acceler- 
ated by radiation pressure, yielding a lower dust abundance. It 
then increases until it reaches a constant value when the ratio 
of v(r) to (v(r) + VdriftM) is constant. A(r) is calculated from 
the adopted gas mass-loss rate M and dust-to-gas mass ratio if/ 
through Eq. (3). 

Considering the motion of grains, and equating the radia- 
tion pressure force with the gas drag force, the drift velocity of 
a grain of size a can be calculated to be (Kwok 1975; Truong- 
Bachetal. 1991) 



Vdrift(a, r) = v K (a, r) [ y/l + x(r) 2 - x(r)] 
with 



(4) 



v K (a,r) = J-r— 
\ M{r)c J 



Q A (a)L A dA 




x(r) 
vr(r) 



with Qa(ci) the extinction efficiency as derived by Justtanont & 
Tielens (1992) for silicates, L A the monochromatic stellar lumi- 
nosity at wavelength A, v T (r) the Maxwellian velocity, T(r) the 
temperature, k the Boltzmann constant, // the mean molecular 
weight of the gas and mu the mass of the hydrogen atom. We 
note that the J94 assumption of Vdrift(a, r) = vk(o, r) holds only 
in the cool outer region where thermal velocities vj are small 
compared to Vdnft- The value of Vdrift(<z, r) should be lower than 
~ 20 km s _1 , at which sputtering of the dust grains starts to be- 
come important (Kwok 1975). When the star is luminous and 
goes through a phase of low mass loss, this value can be higher 
than 20 km s _1 for a certain grain-size a. When this occurs, dust 
grains of size a are destroyed, and hence the contributions of 
the dust grain of size a to the gas-grain collisional heating, the 



sorption cross-sections are in the Rayleigh regime, and hence scale 
with the total dust volume. The gas-grain drag on the other hand scales 
with the geometric cross section. The grain-size distribution also af- 
fects the temperature structure through the gas-grain collisional heat- 
ing. Since smaller grains obtain a smaller drift velocity a steeper grain- 
size distribution causes less efficient heating. On the other hand, pho- 
toelectric yields will be higher for small grains. 



photoelectric heating from dust grains and the heat exchange 
between dust and gas (see Sect. 2.2) is zero. 

Using the above-mentioned equations, the ratio of the drag 
force to the gravitational force T(r) can be written as (J94) 



T(r) = 



I6np, 



3vw_ rr 

: cGM*M(r) J J a[v(r) + v^a, r)] 



2.2. The thermal balance equation of the gas 

From the first law of thermodynamics expressing the conserva- 
tion of energy of the gas, the perfect gas law and the equation 
of mass conservation, the kinetic temperature of the gas in the 
CSE is governed by the relation 



1 dT(r) 
TV) dr 
2 

+ -■ 



3r I 2 



(6) 



H(r) - C(r) 



3 kT{r)v{r)n{U 2 ){r)[Mr) + 1 + / He (r)(/ H (r) + 2)] 
where e(r) = d In v(r)/d In r the logarithmic velocity gradient, 
H(r) the total heating rate per unit volume, C(r) the total cool- 
ing rate per unit volume, fn(r) the number fraction of atomic to 
molecular hydrogen i.e. n(H)(r)/n(H2)(r), and fn e (r) the helium 
abundance i.e. rtHeW/«HW = n(He)(r)/(n(H)(r) + 2n(H 2 )(r)). 
The first term on the right hand side of Eq. (6) represents the 
cooling due to adiabatic expansion in case of constant mass 
loss. The second term represents the balance of the heating and 
collision-driven radiative cooling processes. The denominator 
in this expression takes the contribution of the different colli- 
sion partners H, H2, and He into account. 

In order to calculate the gas kinetic temperature as a func- 
tion of radius, all the major heating and cooling sources must 
be taken into account. The main processes are discussed below. 
Newly implemented energy sources w.r.t. J94 are photoelec- 
tric heating from dust grains, heat exchange between dust and 
gas, and heating by cosmic rays. All the other heating/cooling 
terms underwent updates mainly concerning the used inelas- 
tic collision rates. As demonstrated by J94, cooling by vibra- 
tional excitation of H2O molecules is unimportant, since the 
critical density of H2O is too low in the envelope to play a 
role in the energy balance. We therefore neglect this effect. We 
note that we have not yet accounted for fine-structure line cool- 
ing. This process may be important for supergiants having a 
strong chromosphere (e.g. a Ori, Rodgers & Glassgold 1991). 
A proper calculation of the number densities of the main atoms 
and molecules requires the implementation of a chemical reac- 
tion network programme, which is planned for the future. 

2.2.1. Gas-grain collisional heating 

Both gaseous and dusty material in the CSE are accelerated 
away from the central star by radiation pressure. This radia- 
tion pressure mainly acts on the dust grains, and momentum is 
transferred to the gas by collisions between gas molecules and 
dust grains (GS76). These collisions are the main heat source 
for the gas. The viscous heat input per unit volume can be writ- 
ten as 

1 , 

H gg (r) = n d (a, r)(r d (d)\ A ^(a, r) X ^p(r)v drift (a, r) , (7) 
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with cr d (a) the geometrical cross section of a dust grain of size 
a, being na 2 . Taking the grain distribution into account, one 
obtains 



ffgg(r) = -xAMxm H x n 2 (H 2 )(r)x(/ H (f) + 2) 2 



x(l+4/ He (r))x 



r)da , 



(8) 



with m H the mass of the hydrogen atom. Note that viscous heat- 
ing depends strongly on the value of Vdrift(a, r). 

2.2.2. Photoelectric heating from dust grains 

Simple physical arguments demonstrate that the impact of far 
ultraviolet (UV) radiation (with wavelengths between 912 and 
1200 A) removes electrons from dust grains in such an efficient 
way that about 4 % of radiative energy can be converted into ki- 
netic energy, heating the gas (Tielens & Hollenbach 1985). The 
calculation of the heating rate of the gas by this process is com- 
plicated by the fact that the electron ejection rate and the photo- 
electron energy depend on the charge of the dust grains, which 
is a function of temperature T and electron density n e in the gas. 
For a single size of grain, this can be incorporated into a semi- 
analytical formalism (de Jong 1977). Dealing with a grain size 
distribution complicates the computations, since essentially all 
grains will attain many different ionisation stages under typi- 
cal conditions. Such single and more highly charged grains all 
contribute to the photoelectric heating (Bakes & Tielens 1994). 

Computations for the photoelectric ejection of electrons 
from a dust grain size distribution have only been performed in 
the carbon-rich case. Bakes & Tielens (1994) demonstrated that 
the photoelectric heating rate increases rapidly with decreasing 
grain size, and "classical" grains (a « 0.1 yum) contribute neg- 
ligibly to the total heating rate. Approximately half of the total 
heating rate is due to species with a < 15 A. The other half is 
contributed by larger species up to a < 100 A. This decrease in 
photoelectric efficiency with increasing grain size results partly 
from the decrease in photoelectric yield, i.e. the photoemission 
yield of small particles is expected to be significantly larger 
than that of bulk materials, since for bulk materials the photo- 
electron may lose all its excess energy in collisions with other 
atoms before it reaches the surface. Furthermore, it is aided by 
the increase in average grain charge with increasing grain size, 
which is a direct result of the fact that the UV absorption rate 
increases with grain volume, while the recombination rate de- 
pends only weakly on a. Hence, the lower limit of the grain size 
distribution has a considerable influence on the total photoelec- 
tric heating rate, while the upper limit has only little effect. 

The physical nature of the circumstellar dust grains around 
O-rich stars is far less understood and the coupling between 
the far-UV radiation field and the gas is unclear. To approxi- 
mate the photoelectric heating efficiency, we therefore used the 
analytical expression as obtained by Bakes & Tielens (1994). 
Since this expression was obtained for a grain-size distribution 
with a lower limit of 3 A, we scaled their result with a factor 
0.2 to account for our a m {„ of 50 A(see their Fig. 8) 



#pe(r) = 0.2 x 10~ 24 x 



3 x 10~ 2 



xn(H 2 Xr)(/ H W + 2) (9) 

in erg s _1 ctrT 3 , with Go = 1 the standard value for the intensity 
of the incident far-UV field in units of the Habing interstellar 
radiation field. This equation is obtained for a grain-size dis- 
tribution n d (a,r)da = Aa~ 3 - 5 rmda, with A being 10~ 2516 . The 
electron density n e (r) is calculated from the reaction CO — > 
O + C— >0+C + +e~, which is known to be the main provider 
of electrons in the CSE. As demonstrated by Mamon et al. 
(1988) in their Fig. 6, the abundance of neutral C in the CSE is 
always less than the CO or C + abundance. 

Eq. (9) still has to be scaled with a factor Go x exp(-r uv ), 
where t uv represents the optical depth of the shell against UV 
radiation as measured from the outer radius, i.e. 



#p e (r) = # pe (r) x G x exp(-r uv (r)) . 



(10) 



The optical depth r uv (r) is taken to be 1.8A„(r) (Hollenbach & 
Tielens 1999). In case of an interstellar silicate-to-gas ratio the 
extinction A v — 1 .6 X 10~ 22 x Nh, with Nh the column number 
density of hydrogen. The constant is scaled to the appropriate 
dust-to-gas mass ratio of the shell, i.e. A v = 1.6 x 10~ 22 /0.01 X 

2.2.3. Heat exchange between dust and gas 

Since gas and dust have a different temperature, heat can be 
exchanged between the two species. This process can heat/cool 
the gas considerably, but has almost no effect on the dust tem- 
perature. The heating/cooling rate per unit volume in units of 
erg s _1 ctrT 3 may be written as (Burke & Hollenbach 1983) 



n H (r)n d (a,r)cr d (a) J— — a T (r) 



nniH 



(11) 



x2k[T d (r) - T(r)]da 



= Ank 



8fc 



xA(r)xn 2 (H 2 )(r)(/ H (r) + 2) 2 X 



5F(r) V?V) [ W - T(r)] (a^ 5 - a^) , 



(12) 



1 + 2X 10" 4 G ^JW)/n e (r) 



with T d (r) the dust temperature and a T (r) the average accom- 
modation coefficient (Burke & Hollenbach 1983; Groenewegen 
1994) 

aF(r) = 0.35 exp (- y/(T d (r) + T(r))/500) + 0.1 . (13) 

Since we do not solve the radiative transfer in the dust contin- 
uum, the dust temperature as function of the grain-size T d (a, r) 
is unknown. Instead, the dust temperature is assumed to vary 
as (Olofsson in Habing & Olofsson 2003) 

T d (r) = T*(RJ2r) 2/(4+s \ (14) 

with T* and R* the temperature and radius of the central star. 
Observational data suggest that 5*1 (Olofsson in Habing & 
Olofsson 2003). 



2.2.4. Heating by cosmic rays 

The heating rate per unit volume by cosmic rays is calculated 
by Goldsmith & Langer (1978), and given by 

H a (r) = 6.4 x 10- 28 x «(H 2 )(r)(l + / H (r)/2)(1 + 4/ He (r)) (15) 
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in units of erg s 1 cm 3 . The uncertainty in the numerical factor 
is about a factor of two. 

2.2.5. Cooling by rotational excitation of H 2 

Because the typical temperatures in the CSE around evolved 
stars are low, radiative cooling mechanisms involving the colli- 
sional excitation of the rotational levels of abundant molecules 
form the major sinks of energy. H2O (this section) and CO 
(next section) with their high dissociation energy and substan- 
tial dipole moment are very important in this respect. 

It is beyond the scope of this paper to calculate the pop- 
ulation of the many H 2 rotational levels in a self-consistent 
manner. Instead, we will follow GS76 in approximating the 
structure of a H2O molecule by a three-level system: two ro- 
tational levels in the ground electronic vibrational state v = 
at energy ~ kT and one rotational level at v = 1 . Assuming that 
all of the pure rotational transitions in the ground vibrational 
state have the same excitation temperature r exc , the heat loss 
rate per unit volume per unit time is given by the difference 
between the collisional rate for excitation and de-excitation 



C H2 o,roi(r) = n(H 2 0)(r) x hv 2 ,i(r) x |n(H)(r)<crv>(H-H 2 o)+ 
n(H2)(r)(o-v) ( H 2 -H 2 o) + «(He)(r)((rv>(He-H 2 o)] x 
x[ ex V (-hv v (r)/kT(r)) - exp(-hv 2A (r)/kT ex{: (r))] , (16) 

with V2,i(r) the transition frequency ( = 2.6 x lO n T^ c (r) s _1 ), 
and (o~v)(r) the inelastic collision rate constant. Quite often the 
He-H20 inelastic collisional rate constant is assumed to be a 
factor V2 lower than the H2-H2O rate constant due to the dif- 
ference in mass (e.g. Groenewegen 1994). However, as demon- 
strated by Phillips et al. (1996), excitation by para-H2 is not 
too different from excitation by He atoms at excitation temper- 
atures from 20 to 140 K, with most rates being within a factor 
1-3 larger, but excitation by ortho-H2 is significantly differ- 
ent, with some rates an order of magnitude larger than rates for 
excitation by He atoms. From Tables 6 and 7 in Phillips et al. 
(1996), we use as mean values for the ratio of the inelastic col- 
lision rates 



< o "v) ( H 2 (j=0)-H 2 Op a r a ) 



<0"V> (H c-H 2 O para ) 
<0"V>H 2 (J=1)-H 2 0p ala 



<C r v)(He-H 2 O para ) 



1.7 



*7.7 



<0"V>H 2 (j=O)-H 2 O o , 



<0T>(He-H 2 O ortho ) 
<CrV>H 2 (j=l)-H 2 OBho 



<crv> (H< 



e-H 2 or , ho ) 



1.6 



*7.8. 



Theoretical rate constants among the lowest 45 para and 45 
ortho rotational levels of water in collisions with He atoms have 
been calculated by Green et al. (1993). From their Tables 2 and 
3, we derive 



<crv> ( He-H 2 o para ) ~ 0.21 x y/T(r) x 10" 11 cm 3 sec -1 
<<™>(He-H 2 cw ~ 0.13 x x 10- 11 cm 3 sec" 1 . 

The H-H 2 inelastic rate constant is assumed to be a factor 
1.16 larger than the He-H 2 constant (this value takes both the 
smaller mass and the smaller cross section of H into account). 
With an ortho-to-para ratio of 3:1 for both H2 and H2O, we 
obtain as cooling rate for rotational excitation of H2O 



Cu 2 oUr) = «CH 2 )(r) n(K 2 0)(r) hv 2il (r)[0.2\ x 10" n x 
Vrw] x [0.83/hW + 0.715/ He (r){/ H (r) + 2) + 4.5] x 
x[exp{-hv v (r)/kT(r)} - ex V {-hv 2il (r)/kT exc (r)}] . (17) 

The excitation temperature r exc (r) can be calculated from 
the rate equations, using an escape probability formalism to 
describe the radiative transfer in the line; see Eq. (1 1) in J94. 

The water molecules in the CSE are photodissoci- 
ated by interstellar ultraviolet photons with A « 1650 A. 
In particular, the photochemical chain is given by 
H 2 ->OH + H^O + H + H (GS76; Huggins & Glassgold 
1982; Netzer & Knapp 1987). This photodissociation of H 2 
is a significant source of OH molecules in the outer envelope. 
Since no functional form exists for the distribution of water 
in the CSE — and the implementation of a chemical reaction 
network programme is planned for the future — we assume a 
similar description for the H2O abundance profile as derived 
by Mamon et al. (1988) for CO (see next section, Eq. (22)). 
In this functional form, we have to provide the position r^ ' 
where the abundance of water relative to H2 has decreased by 
a factor of two relative to the value at the stellar surface. This 
value is set by (Groenewegen 1994) 



,(H 2 0) _ 



'1/2 



= 35 



/ — \ - 7 
M/10" 5 Moyr -1 ) 



x(v(O/15kms -1 r a4 10 15 cm, 



(18) 



where (in case of variable mass loss) M is taken as the mean of 
M(r). This expression is derived from the results of Netzer & 
Knapp (1987), who determined the OH peak number densities 
in the CSEs around evolved stars. The thus assumed H2O dis- 
tribution is similar to the results of Nejad & Millar (1988, see 
their Fig. 1). 

2.2.6. Cooling by rotational excitation of CO 

Justtanont et al. (1994) included CO cooling in a similar way 
as H2O rotational cooling. The H-CO inelastic rate constant is 
assumed to be a factor 1.16 larger than the He-CO inelastic 
rate constant due to the difference in mass (velocity) and cross 
section. Using (Schinke et al. 1985) 



<0-v> ( co-H 2 ) * 4 x !0~ 12 x cmV 1 , 

and 

<CTV> ( C0-He) ~ 4.5 X 10~ 13 X CmV 1 , 

(McKee et al. 1982), we obtain 

C CO ,rotW = «(H 2 )(r) x «(CO)(r) x hv 2 ,i{r) 
X [4.5 x 10~ 13 x ^Jf{r)\ 

X [1.16/hW + /HeWl/HW + 2} + 10] 

X [cxp{-hv 2 ,i(r)/kT(r)} - exp{-hv 2 ,i(r)/kT x (r)}] . 



(19) 



(20) 



(21) 
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with v 2 ,i(r) = 4.89 x 10 10 r A 05 (r) s"\ and T x (r) the excitation 
temperature for CO molecules defined in a similar manner as 
H 2 (see previous section). 

An important quantity is the CO abundance in the photo- 
sphere, which is determined by the carbon and oxygen supply. 
During the first dredge-up on the red giant branch the carbon 
abundance is depleted to roughly two-thirds of the main se- 
quence value while the oxygen abundance remains almost un- 
changed. On the AGB, carbon may be added to the envelope 
due to the third dredge-up. Typical CO and H 2 abundances 
at the base of the CSE are f C o(R*) = n(CO)(tf*)/n(H 2 (#*)) = 
5.3xlO- 4 and/H 2 oOR*) = n(H 2 0)(/?*)/n(H 2 (fl*)) = 1.2 xlfr 3 
when the star arrives on the AGB. If the C/O ratio in an AGB 
star has increased to 0.7, the abundances changes to 1.2 x 10~ 3 
and 5.1 x 10 4 for CO and H 2 respectively (Groenewegen 
1994). Knapp & Morris (1985) arrived at a value of f C0 (R*) = 
3 x 10~ 4 for O-rich giants, 6 x 10~ 4 for S-type stars, and 8 x 10~ 4 
for C-rich giants; Zuckerman & Dyck (1986) used 5 x 10~ 4 for 
O-rich and 1 X 10~ 3 for S- and C-rich giants. This uncertainty in 
fco translates in an uncertainty on the derived mass-loss rate. 
In our code the carbon and oxygen abundances at the stellar 
radius are free input parameters, and can be adjusted ac- 
cording to the evolutionary state of the target under study. If 
no information can be found on the carbon and oxygen abun- 
dances in the outer atmosphere of the target, we opt to use the 
cosmic abundance values of Anders & Grevesse (1989). 

CO will be dissociated by the interstellar UV field, and its 
abundance will decrease outward. Mamon et al. (1988) have 
investigated this effect on the CO abundance distribution. They 
obtained that the CO abundance f(co)( r )> relative to the value 
at the photosphere /(co)(R*)> can be approximated by 



/(CO)W = / (C 0)(R*)exp[-ln2(r/^ 2 0) ) ff ] , 



(22) 



where r?^ denotes the position where the number density has 



1/2 



decreased by a factor of two. The values for rfj^ and a de- 
pend on the gas velocity and mass loss rate and are tabulated 
by Mamon et al. (1988). 

2.2.7. Cooling by vibrational excitation of H 2 

Cooling by vibrational excitation of H 2 molecules is simpler to 
calculate than for H 2 0. The reasons are: (1) The vibrational 
energy level spacing is large (~ 0.6 eV). Hence only colli- 
sional excitation of the first excited vibrational state needs to 
be considered. (2) Since H 2 is dipole inactive, resulting in small 
Einstein A-values, the lines are optically thin even in the dense 
part of the CSE. (3) H 2 is in LTE through a large part of the 
envelope (GS76, J94). The description of the heat loss per unit 
volume is taken from GS76, and can be summarised as 



Ch 2 = A h0 hvi. m(r), 



(23) 



where A^o ~ 3 x 10~ 7 s~' is the spontaneous emission rate 
from the first vibrational state, hv\$ as 0.6 eV is the energy of 
the emitted photon, and tt\ - n H2 (v = 1) the number density 
of vibrationally excited H 2 molecules. The number density n\ 



of vibrationally excited H 2 by collisions with H 2 and H atoms, 
can be written as 

f ^ ru w n a„(r) X exp(-/tvi, /W(r)) 
ni(r) = n(H 2 )(r) — , , r , . ______ , „ . ; — (24) 



a n (r) x [1 + exp(-hv hQ /kT(r))] + A uo 



with 



a n (r) = n(H)(r)<cr*Vrt> (H -H 2 ) + n(H 2 )(r)<o-*v, A > (H2 -H 2 ) (25) 

where (<x*v,/,)(h-h 2 ) and (cr*v t h) (h 2 -h 2 ) are the collisional de- 
excitation rate constants. 

Theoretical values for <o"*v,/,)(h-h 2 ) are approximated by 
(Hollenbach & McKee 1979) 

<cr*v, ft ) (H -H 2 ) = 1.0 X 10- 12 7\r) 1/2 x exp [-(1000/7\r))] (26) 
cm 3 s _1 , and for <0-*v,/,} ( h 2 -h 2 ) (Hollenbach & McKee 1989) 
<ct*v, a ) ( h 2 -h 2 ) = 1.4xl0- 12 7Xr) 1/2 



xexp{-[18100/(r(r)+ 1200)]} 



(27) 



3 -1 

cm s . 



2.3. Treatment of the boundary conditions 

The differential equations governing the temperature and ve- 
locity structure (Eqs. (2) and (6)) are solved simultaneously. 
For the central star, we use a blackbody with a user-defined 
temperature T*. As inner boundary for the velocity structure, 
we assume that the flow velocity of the gas is equal to the lo- 
cal sound velocity, v s , at the radius of the CSE, Ri nner , where the 
dust condenses. Since many questions still exist on the complex 
structure between R+ and Ri nner we assume that the temperature 
in this regime is described by a power law 



T(r) = T^ 



(28) 



with £ « 0.5, which holds when both emission and absorption 
are optically thin. The velocity law can be approximated by the 
classical /3-law 



v W ^(l-^f 
with 



Ro = R*U 



VP 



(29) 



(30) 



with vo the velocity at the photosphere, which can be calcu- 
lated directly from v s , and (3 being 1/2, a typical value in case 
of winds of cool stars (Schutte & Tielens 1989). Since we 
are studying low excitation rotational CO line profiles being 
formed quite far in the envelope, the exact value of e and f3 
does not influence the calculated line profiles. The outer radius 
is determined by the CO abundance, and is set at the radius 
where the CO abundance fco( r ) drops to 1 % of its value at 
the photosphere. Adopting a value for the gas mass-loss rate, 
the dust-to-gas ratio is varied until the observed terminal ve- 
locity is obtained. Since F(r) oc if/ ■ Q A (a) (Eq. 5), the derived 
dust-to-gas ratio if/ should be interpreted in relation to the used 



L. Decin et al.: Probing the mass-loss history of AGB and red supergiant stars from CO rotational line profiles 



7 



tabulated extinction efficiencies Q A (a) derived by Justtanont & 
Tielens (1992) for silicate grains. 

The sudden onset of some heating or cooling terms, yields 
a change in going from a non-stiff to a stiff set of differen- 
tial equations (or vice versa), which a standard Runge-Kutta 
method is not able to handle. We therefore chose to use the 
DLSODA package as provided by ODEPACK 2 . Going from 
the inner boundary to the outer boundary, this solver auto- 
matically selects between non-stiff and stiff methods. It uses 
Adams methods (predictor-corrector) in the non-stiff case, and 
Backward Differentiation Formula methods (the Gear meth- 
ods) in the stiff case. We require an absolute accuracy of 1 K in 
temperature and lOOcms -1 in velocity, i.e. the estimated local 
error in the parameter is less than the above-mentioned accu- 
racy. 

2.4. Variable mass loss 

It is beyond the scope of this project to implement a full hy- 
drodynamical code to predict variabilities in the wind struc- 
ture. We simulate the effect of various types of variable mass 
loss, as e.g. the superwind phase or "quasi-periodic variations" 
discussed in Sect. 1, in an empirical manner by adapting the 
mass-loss rate as a function of radius, i.e. M(r). 

In order to incorporate a variable mass loss, we calculate for 
each radial point r the expected temperature and velocity as if 
resulting from a constant mass outflow with the corresponding 
mass loss rate. In practise this means that we solve Eqs. (2) and 
(6) from the inner edge Ri nner up to r with a fixed M chosen 
to correspond to M(r). In this way we assume that the main 
influence on the temperature is set by the expansion of the flow 
and the local heating and cooling terms. In this approach the 
influence of the neighbouring shells, which may correspond to 
different mass-loss rates, is neglected. 

If the mass-loss behaviour is changed, the drift velocities 
will be affected. This in turn causes a change in the flow ve- 
locity v(r) of the gas. This feedback is self-consistently ac- 
counted for in the calculations of the heating and cooling rates. 
However, we do not take these modified velocities into account 
when solving the radiative transfer equation, because our the- 
oretical code solves the radiative transfer using the comoving 
frame (CMF) method as developed by Mihalas et al. (1975, see 
next section) which demands that the velocity should increase 
continuously with radius. In case of variable mass loss, we use 
the velocity structure as inferred from the case with constant M. 
For instance, in the test case of GX Mon described in Sect. 2.5, 
the velocity structure as in the absence of a superwind is used. 

Several options are built in the code to define M(r), be- 
ing constant, a heavy-side function, a block function, a picket- 
fence function or analogous functions with smoother deriva- 
tives. 

2.5. Discussion 

To check the new input of heating and cooling rates, and new 
numerical routines we simulated the results as obtained by J94 

2 http://www.netlib.org/odepack/ 



Table 1. Parameters of GX Mon as derived by J94. 



parameter 


value 


T* 


2500 K 


R* 


4x 10 13 cm 


M, 


1M 


Rinner 


6.25 R, 


M 


7.2 x 10- 6 M o yr' 


v stoch 


1 km s -1 


distance 


740 pc 




19kms' 



for the Mira-type star GX Mon. Parameters as derived by J94 
are listed in Table 1 . Our final temperature and velocity struc- 
ture are displayed in Fig. 1 in the lower right and upper right 
panel respectively, and should be compared with Fig. 1 in J94. 

Inspecting the gas and drift velocity structure shows quite 
a difference between both results. The reason for this differ- 
ence is a mistake in the implementation of the results of the 
ordinary differential equation (ODE) operator (being a Runge- 
Kutta method with adaptive stepsize control) in the version of 
J94. J94 erroneously assigned the computed parameters being 
the velocity v(rl), temperature T(rl) and the respective deriva- 
tives dv(rl)/dr and dT(rl)/dr at radius point rl also to the next 
radius point r2 in the grid, which input was then used to calcu- 
late the parameters at the newly derived radius point r2', com- 
puted using adaptive stepsize control. This inaccuracy yielded 
in case of GX Mon a terminal velocity being ~ 8 km s _1 too 
low and a terminal drift velocity being ~ 5 km s _1 too low. 
Moreover, our improved treatment of the drift velocity (Eq. 4, 
Sect. 2.1) lowers for GX Mon the terminal drift velocity by 
0.5kms _1 and the terminal flow velocity by 2kms _1 com- 
pared to the approximation of Vdrift(a, r) = v K (a, r) (J94). 

That there is some difference between the derived temper- 
ature structures does not come as a surprise, since new heating 
terms are included and some of the ones already accounted for 
underwent major updates. 

- Heat exchange between dust and gas and cosmic ray heat- 
ing have been included, but they only have a minor contri- 
bution to the total heating rate. 

- The gas-grain collisional heating rate as given by J94 is 
somewhat too low due to a numerical mistake in the imple- 
mentation of the fractional abundance of atomic hydrogen 
/h yielding too low a scaling factor A. Since this frictional 
heating is the main heating source over most of the enve- 
lope regime, our calculated T(r) is in general higher than in 
J94. 

- The main differences in the cooling rates come from the up- 
dates of the inelastic collisional rates involved in the cool- 
ing by rotational excitation of H2O and CO, being lower 
than in J94. 

As an illustration of the effects of variations in the mass- 
loss rate, we simulate a superwind phase with an enhanced 
mass loss of a factor of ten till 200 R* ( = 8x 10 15 cm) compared 
to the parameters of the base-line model of GX Mon given in 
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Figure 1. The velocity structure and gas kinetic temperature from calculating the energy sources and sinks in the CSE of GX Mon. The start of 
the dusty envelope, Rj nlrer , is indicated by the dotted line. The drift velocity is displayed for the maximum grain-size, a = 0.25 /im. Note that for 
grains with radius smaller than the wavelength, Q oc a; hence vj ri f t oc a 05 . First panel: Cooling rates; second panel: heating rates; third panel: 
velocity structure; fourth panel: gas kinetic temperature. 



Table 1 . We do not attempt to better explain the observed line 
profiles of GX Mon. In Sect. 4 we will present the results of a 
parameter study focusing on explaining the observations of VY 
CMa. 

Fig. 2 displays the temperature, velocity, heating and cool- 
ing terms for this simulation. As explained in Sect. 2.4, the ve- 
locity structure used in the radiative transfer calculations is not 
coupled to a change in mass-loss rate, but is kept as inferred 
from a case with constant M. In the panel giving the velocity 
structure, the black lines represent the gas and drift velocity in 
case M- 7.2 X 10~ 5 M yr~' would have been used to calculate 
the velocity structure used in the radiative transfer calculations, 
the blue (grey) line is for M- 7.2 x 10~ 6 M Q yr _1 . An increase 

. -1/2 

in M yields (1) a decrease of the drift velocity (vadft <* M , 
Eq. (4)) since the coupling between dust and gas increases, and 
(2) a higher terminal velocity, since there is more dust to inter- 
cept the radiation and hence to accelerate the gas by transfer of 
momentum. Which constant M value to choose in this case can 
be debated, but as illustrated in Fig. 2, the difference between 
both v(r)-structures is small. A more correct gas velocity pro- 
file should display a decrease of the flow velocity from 200 R* 
( = 8 x 10 15 cm)on. 



The change in viscous heating is not dominated by the de- 
crease in drift velocity, but by the increase in dust and H2 par- 
ticles, yielding a higher gas-grain collisional heating rate out 
to 200 R*. Other energy rates are increased too, with the total 
influence being quite complex. In general, the temperature is 
higher till 200 R* w.r.t. the constant mass-loss case displayed 
in Fig. 1. The effects of this superwind phase on the line pro- 
files are discussed in Sect. 3.3. 

3. Line radiative transfer 

Late type giants and supergiants have a ratio of stochastic to ter- 
minal velocity that may be in the order of magnitude of ~ 0.5 
or even 1 (Che et al. 1983) yielding rather broad scattering 
zones and thus a highly non-local radiative transfer. Therefore 
the condition for using the Sobolev approximation (Sobolev 
1947; Castor 1970) is not fulfilled anymore. Schonberg (1985) 
demonstrated that non-local scattering effects cause severe er- 
rors in the source function when using the Sobolev approxima- 
tion (see also J94). 

Moreover, as demonstrated by Schonberg (1988), part of 
the line profile of the low-excitation CO transitions is formed 
under non-LTE circumstances. The decoupling radius (i.e. that 
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Figure 2. Same as Fig. 1, but now with a superwind phase, i.e. with a mass-loss rate being a factor of 10 higher up to 200 R* (indicated by the 
dashed line). The velocity structure displayed in Fig. 1 is also shown in grey in the third panel. More details are given in the text. 



distance from the star where the excitation temperature in the 
lines is by ten percent smaller than the kinetic temperature) is 
a function of the specific line. The higher the excitation level, 
the smaller the decoupling radius. Clearly, the decoupling ra- 
dius decreases with decreasing mass-loss rate. For a model 
with parameters almost identical to the ones listed in Table 1, 
Schonberg (1988) demonstrated that the excitation of the rota- 
tional 12 CO lines in the vibrational ground state is dominated 
by collisions in the mass-loss range M k, 3 x 1CT 5 M yr _1 . 

We therefore opt to calculate the CO rotational line profiles 
using a non-LTE radiative transfer code based on the multilevel 
approximate Newton-Raphson (ANR) operator as derived by 
Schonberg & Hempe (1986). The ANR operator has been de- 
rived from the comoving frame equation of radiation transfer 
(Mihalas et al. 1975). The method is similar to the partial lin- 
earision approach. The ANR operator is a diagonal (therefore 
local) operator, and can be derived with only a slight increase 
in computational effort over that required to find the formal so- 
lution. When solving for multilevel line formation in expand- 
ing envelopes, Schonberg & Hempe (1986) have shown that 
convergence was reached within 10-20 iterations for winds 
around cool giants. We note that although the method is proven 
to be stable, Hillier (1990) demonstrated that at very high opti- 
cal depths the diagonal ANR operator may give rise to damped, 
depth-coupled oscillations, which may be due to the inability of 



the diagonal operator to correctly handle the diffusion approxi- 
mation. He has shown that in these cases tridiagonal (or penta- 
diagonal) Newton-Raphson operators have a superior conver- 
gence rate. 

A description of the radiative transfer equations can be 
found in Schonberg (1988). The calculations include excitation 
by collisions with H2 as well as scattering of infrared and mi- 
crowave radiation emergent from warm dust particles, from the 
stellar photosphere and the cosmic background radiation. The 
implementation of the dust opacity was modified to include the 
grain opacity of silicates as computed in Justtanont & Tielens 
(1992). 

After solving the radiative transfer, line profiles for each 
transition are calculated by ray-tracing using formal integration 
(Schonberg 1988). Writing the calculated flux as F v , the main 
beam temperature T M b is calculated by (Schonberg 1988) 

r MB = ^v-^, (31) 

with Djel the diameter of the antenna. Note that the line pro- 
file will not be symmetric since the opacity in the blue wing is 
produced at slightly larger r than in the red wing. This kind of 
asymmetry in double-peaked (optically thin) line profiles will 
be stronger for resolved envelopes, and provides a useful di- 
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Figure 3. Comparison of the benchmark case la (left) and lb (right) as defined in van Zadelhoff et al. (2002) for different spatial resolutions. 
Upper panel: number density n 2 in the upper level of the Active 2-level molecule as a function of radial distance from the star. The full line 
represents the mean of the number densities of all codes discussed by van Zadelhoff et al. (2002). Lower panel: ratio between their and our result. 
For increasing resolution the solutions converge. Dotted line: AR/R ss 0.07; dashed line: AR/R ss 0.02; dashed-dotted line AR/R « 0.005. 



agnostic for obtaining the stochastic velocity (see Schonberg 
1988). 

3.1. The molecular model 

The rate coefficients for collisions of CO with H2 at kinetic 
temperatures from 10 to 4000 K are taken from Larsson et al. 
(2002). In this work, rotational quantum numbers up to /„ = 40 
were taken into account. Data on the vibrational-rotational 
v = 1 - and the rotational transitions are from the HITRAN 
database (Rothman et al. 1987). A total of 62 levels and 120 
transitions are treated, of which 31 levels are in the ground 
state of 12 C 16 0, and 31 levels in the first vibrational state. 
Comparison with the CO line list of Goorvitch & Chackerian 
(1994) shows that the used Einstein A-coefficients are in gen- 
eral good agreement; they are however systematically lower 
than the ones of Goorvitch & Chackerian (1994), with a mean 
difference of a factor of 1.017, a maximum of a factor 1.035, 
and a minimum of a factor 1.005. 

3.2. Benchmarking 

To test the accuracy, the new code has been benchmarked 
against existing radiative transfer codes. We used the study per- 
formed by van Zadelhoff et al. (2002), who compared a num- 
ber of independent computer programs for radiative transfer in 
molecular rotational lines. None of the codes discussed in this 
comparison solves the radiative transfer in the comoving frame. 
While in the rest (inertial) frame the precision of a numerical 
code can be increased by using a second-order (or higher) inte- 
gration of the transfer equation, in the comoving frame formal- 
ism it is written as two first-order differential equations yield- 
ing only first-order accuracy (Mihalas et al. 1975). In the case 
study of problems la/lb [having a simple power law density, 
constant temperature and a Active 2-level molecule], the use of 
a first-order integration introduces errors of the order of 12 % 



(van Zadelhoff et al. 2002). Hence, for a proper benchmarking 
of the code, the spatial resolution in the CMF should be in- 
creased to 'mimic' a higher-order solution, i.e. to improve the 
accuracy. This is illustrated in Fig. 3, for case la and lb. It is 
clear that higher spatial resolution results in a better agreement 
with the other (non-CMF) codes, reaching a 1.2 % level of ac- 
curacy for AR/R * 0.005. 

Case 2 of van Zadelhoff et al. (2002), describing an inside- 
out collapsing envelope observed in rotational transitions of 
HCO + , has not been simulated since the CMF method is de- 
signed for a continuous increasing velocity, implying a blue- 
wing condition stating that the highest frequency point in the 
local profile can intercept only continuum radiation from other 
points in the atmosphere. 

3.3. Effect of variable mass loss on the predicted line 
profiles of CO 

For the specific case of GX Mon and the simulation of a super- 
wind phase with the same parameters as in Sect. 2.4 the cal- 
culated low excitation CO line profiles are displayed in Fig. 4. 
The adopted beam sizes are those for the JCMT, and are in- 
dicated in the upper corner of each panel. Number densities, 
source functions S v and absorption coefficient k v as a function 
of radius are plotted in Figs. 5 and 6. 

A superwind phase increases the number densities of all the 
energy levels from R mna to 200 R*, where the high mass-loss 
steps. 

The sudden rise in temperature at the start of dust conden- 
sation gives a strong peak in the source function S v at Rmner- 
The higher temperature in the superwind phase, results in (1) 
higher values of the source function up to 200 R+, and lower 
values thereafter, and (2) a strong increase in the opacity up to 
200 R+ due to the higher number densities till the sudden de- 
crease in mass loss rate introduces a sharp decrease in opacity. 
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Table 2. Parameters of VY CMa. We have assumed cosmic carbon 
and oxygen abundances (Anders & Grevesse 1989). 



20 -20 20 
Velocity [km s" 1 ] 

Figure 4. CO rotational line profiles (a) full line: for GX Mon using 
T(r) and v(r) of Fig. 1, (b) dotted line: in case of a variable mass-loss 
rate as displayed in Fig. 2. Note the different ordinate axes between 
upper and lower panels. 



The figures displaying I{p) * p 3 are instructive plots to il- 
lustrate where the lines mostly originate (see also Kemper et al. 
2003). In case of a constant mass-loss rate, there is only one 
peak in I(p) * p 3 , being closer to the star for higher excita- 
tion lines; for the superwind case, we see a second peak arising 
around 200 R*, the end of the enhanced mass loss. A direct 
result from the enhancement in number densities at distances 
< 200 R* is an increase in integrated intensity for the super- 
wind case, especially for the higher excitational rotational lines 
which are formed closer in. Moreover, the sudden change in 
opacity (number density) at 200 R* gives rise to 'bumps' in 
the line profiles, here with (projected) velocities around 10 and 
15 km s _1 . 

Depending on the stellar input parameters, the abruptness 
of the variation in M, and the place where (or time when) this 
change has taken place different kind of 'composite' profiles 
(i.e. displaying "dips and bumps" in the line profile) are seen. 
This kind of composite profiles is also seen in observational 
data, with a Ori and VY CMa being examples in Kemper et al. 
(2003), and other targets in Winters et al. (2003). Note that the 
relative strength of the observed profiles of some stars, where 
this kind of 'composite' profile is not visible, is also indicative 
for a variability in the mass-loss rate (see the discussion on WX 
Psc in Kemper et al. 2003). 

The complex interplay between stellar and mass-loss pa- 
rameters and some general trends linking the mass-loss rate to 
the observed line intensity will be described in a forthcoming 
paper. 



Parameter 


This work 


Smith et al. (2001) 


T* 


2800 K 


3000 K 


R* 


1.6 x 10 14 cm 


2x 10 13 cm 


L* 


3 x 10 5 L 


1.6 x 10 5 L o 


M, 


15 M 


15 M 


D 

-^Mnner 


10 R, 


7R* 


M 


see Fig. 7 


10~ 5 x lO^Moyr 1 till 0.5" 


M 


see Fig. 7 


3 x 10~ 4 M yr 1 from 0.5" up to 7 


Vstoch 


1 km s -1 


1 km s -1 


Voo 


35kms _1 


35kms~' 


distance 


1500 pc 


1500 pc 




0.002 


0.01 



tational line transitions emitted in the CSE around VY CMa. 
In Sect. 4.1, the basic parameters of this supergiant are sum- 
marised. The CO line profiles used in the analysis are described 
in Sect. 4.2. Sections 4.3 and 4.4 are devoted to the modelling 
of the CO profiles. We will demonstrate that it is possible to 
derive the mass-loss history of this target in great detail. Our 
results are compared to other studies in Sect. 4.5, and discussed 
in Sect. 4.6. 

4.1. Basic parameters of VY CMa 

VY CMa is a highly luminous, variable M2/3II supergiant. The 
star is heavily obscured, with only 1 percent of the total lu- 
minosity being detected at optical wavelengths (Le Sidaner & 
Le Bertre 1996). The distance, luminosity, effective tempera- 
ture and stellar radius have been estimated by several groups. 
These studies seem to converge on a distance of 1500pc, a lu- 
minosity between 2 - 5 x 10 5 L , and an effective temperature 
of 2800 K (Harwit et al. 2001; Le Sidaner & Le Bertre 1996; 
Monnier et al. 1999). The current mass of VY CMa is estimated 
to be 15 M (Wittkowski et al. 1998). These values are adopted 
for the stellar parameters of VY CMa; they are summarised in 
Col. 2 in Table 2. 

The OH, H2O, and SiO maser emission indicate ex- 
pansion velocities between ~ 30 and ~ 40kms~ 1 (Reid & 
Dickinson 1976; Snyder & Buhl 1975; Reid & Muhleman 
1978). Modelling of CO rotational lines, OH maser emission 
and (near)-infrared data tracing the dust properties indicate val- 
ues for the mass loss between 2.3 X 10~ 5 and 4 x 10~ 4 M yr~' 
(see Table 3). VY CMa's very high mass-loss rate fuels its 
optically thick circumstellar envelope and creates a rich in- 
frared excess spectrum (see the 2.4 - 200 pm ISO spectrum in 
Harwit et al. 2001): metallic Fe and amorphous silicates dom- 
inate the infrared spectrum, with possible minor contributions 
from crystalline forsterite and crystalline enstatite. 



4. VY CMa 



4.2. Observational data of VY CMa 



To illustrate the diagnostic strength of the rotational CO pro- 
files, we now apply our code to model the observed CO ro- 



The rotational CO (2-1), (3-2), (6-5) and (7-6) line profiles 
have been observed by Kemper et al. (2003) with the JCMT 
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Figure 5. Detailed look at the wind structure of GX Mon for a constant mass loss of 7.2 x 1(T 6 M Q yr~' . Temperature [K], normalised number 
density [cm 4 ], source function [erg s" 1 cm" 2 Hz 1 ster -1 ], and the absorption coefficient k v [cm -1 ] as a function of radius; I(p) * p 3 as function 
of the impact parameter p — using T(r) and v{r) of Fig. 1. 



Table 3. Mass-loss rates for VY CMa as found in literature and scaled to a distance of 1500 pc. 





Observ. data 


MtMeyr- 1 ] 


comments 


Reference 




OH maser emission 


1.2 x 10~ 4 


from OH shell radius + n H = 10 4 cm 3 


Booth etal. (1984) 




CO (2-1) 


1 x 10- 4 


scaling based on analysis of 


Zuckerman & Dyck (1986) 








Knapp& Morris (1985) 




CO(1-0) 


2.3 x 10" 5 


from analytical expression 


Loup etal. (1993) 


% 


spectrosc. data 2-2.5 pm + 
photometr. data till 20 pm 
submm data at 400 pm 
ISI 11 /mi interf. 
spectrosc. data 2.38 - 200 pm 
surface photometry 


2.3 x 10~ 4 

3.1 x 10~ 4 
4 x 10~ 4 
10- 5 - 10- 4 
3 x 10~ 4 


data from Hyland et al. (1972); 
xfi = 0.01 

Mdu St = 7.8xlO- 7 M Q yr- 1 
i// = 0.005 
<fr = 0.01 
till 0.5" 

from 0.5" up to 7" V ' 


Bowers etal. (1983) 

Sopka et al. (1989) 
Danchietal. (1994) 
Harwit etal. (2001) 
Smith etal. (2001) 



and the CO (1-0) line by Nyman et al. (1992) with the SEST. 
The flux calibration accuracy for the JCMT data was estimated 
to be 10% for the CO (2-1) and (3-2) lines, and ~ 30% for 
the CO (6-5) and (7-6) lines since no reliable flux standards 
are available for these higher excitation lines. The CO(1-0) 
line has a main beam temperature of 0.06 K, with an rms noise 
of 0.019 K (Nyman et al. 1992). As can be seen from Fig. 8, 
the profiles show a complex structure. Especially the different 
peaks present in the CO (2-1) and (3-2) profiles are indicative 
for variation(s) in the mass-loss rate. 



4.3. A spherically symmetric envelope? 

There exist indications that the CSE surrounding VY CMa is 
not spherically symmetric. While some authors interpret the 
complex optical and infrared data in terms of an expanding 
disk or a bipolar outflow near the central object (< 0.2"; 
e.g. Herbig 1972; Bowers et al. 1983; Efstathiou & Rowan- 
Robinson 1990), other data do not reveal indications for a disk- 
like or bipolar geometry (e.g. Monnier et al. 1999). Recently, 
"Hubble Space Telescope " in combination with ground-based 
infrared images revealed a complex structure of knots and fil- 
amentary arcs in the asymmetric reflection nebula several arc- 
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Figure 6. Detailed look at the wind structure of a superwind phase for GX Mon. The superwind increases the mass loss by a factor 
the stellar surface up to 200 R*. Temperature [K], normalised number density [cm -3 ], source function S y [ergs -1 cm -2 Hz -1 ster - '], 
absorption coefficient k v [cm -1 ] as a function of radius; I(p) * p 3 as function of the impact parameter p. 



10 from 
and the 



seconds across, around the obscured central star (Smith et al. 
2001). 

Nevertheless, we will assume spherical symmetry in mod- 
elling the molecular lines arising from the CSE. This is jus- 
tified with the following arguments: (1) If our assumption of 
spherical symmetry is incorrect, this is more so for the inner- 
most regions, where there are indications of a disk with extent 
£ 0.2" (~ 30 R*). We are however aiming at a reconstruction of 
the mass-loss history from low-excitation CO lines, which are 
formed further out, and trace a region up to a few thousand stel- 
lar radii. (2) If the inhomogeneities are due to blobs of which 
the dimensions are modest compared to the total extent of the 
wind and which are more or less homogeneously distributed in 
the wind, the use of a spherically symmetric model gives us 
the average of the mass loss over all directions. In that way, 
our spherical model gives us an idea on the general mass-loss 
history of VY CMa. (3) Since the beam profiles of the instru- 
ments used to observe the CO lines are quite large (with a half 
power beam width (HPBW) of 45" in the CO(1-0) line, 19.7" 
in the CO(2-l), 13.2" in the CO(3-2), 8.0" in the CO(6-5) and 
6.0" in the CO(7-6) line) the CO line profiles can not be used 
to trace small-scale structure, but only to estimate the density, 
temperature and velocity structure averaged over all directions 
in the regions where they are formed. 

Moreover, Smith et al. (2001) demonstrated that the ap- 
parent large-scale asymmetry of the nebula in the optical and 
near-IR images results from a combination of high extinction 



and backscattering by dust grains. Thermal-infrared images 
(5 - 10//m) reveal a more symmetric distribution. 

4.4. Modelling the rotational CO line profiles 

To determine the mass-loss history of VY CMa a first initial 
grid of ~ 20000 models, with stellar parameters as specified 
in Col. 2 in Table 2, was run to find the global minima in 
parameter- space using a classical reduced ^-method. In this, 
we treat the integrated strength of the line and the line shape 
as two separate entities to be fitted. This allows to give more 
weight to the line profile which is more reliable than the ab- 
solute flux calibration. The mass-loss parameters of the best 
models were then further refined. 

Parameters characterising the CSE of the best model are 
given in the second column of Table 2. Since the calibration 
uncertainty on the CO(2-l) and (3-2) line is smaller than for 
the other profiles, highest weight is given automatically to these 
lines in the reduced ^ 2 -test. The derived temperature profile, 
velocity structure and mass-loss history are displayed in Fig. 7; 
a comparison between observed and theoretical line profiles is 
shown in Fig. 8. 

Before comparing to results of other studies, we first dis- 
cuss the derived M-profile and the associated uncertainties. The 
M-profile has a value of 8 X 10" 5 M yr -1 till 80 R* (denoted as 
(A) in Fig. 7). Mass-loss rates in the order of ~ lxlO -4 M yr -1 
in these inner regions still result in acceptable line profiles for 
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Figure 7. Upper: Estimated temperature profile, middle: estimated ve- 
locity structure, and bottom: estimated mass-loss rate M(r) for VY 
CMa as a function of radial distance to the star (black line). The 
blue/grey line in the panel displaying the mass-loss history repre- 
sents a model with an extra enhancement in mass loss (M(r) = 
2 x 1(T 5 M Q yr" 1 ) around 300 R* having the same theoretical line pro- 
files as in Fig. 8. 



the CO(6-5) and (7-6) line. Between ~ 80 and ~ 600 R*, 
we determine a low mass-loss rate, being of the order of 
1 - 3 x 10~ 6 M Q yr' (denoted as (B) in Fig. 7). This mass- 
loss value can not be below 5 x 10~ 7 M Q yr -1 , since the total 
dust content at the base of the wind would be too low to drive 
the stellar wind. This period of low mass-loss rate was pre- 
ceded (in time) by a short period (between 80 and lOOyr at 
a expansion velocity of 35 km s~') of very strong mass loss 



Figure 8. CO rotational line profiles of VY 
CMa (Kemper et al. 2003) (grey line) com- 
pared with the model predictions (black 
line) based on the parameters of our best 
model (see Table 2). 

in the order of 3 - 3.5 x lO^Moyr 1 (region (Q). This pe- 
riod is situated around 680 R* (~ 1000 yr ago). One can sim- 
ulate almost analogous line profiles with a somewhat smaller 
M (around 2.5 x 10~ 4 M Q yr -1 ) having lasted somewhat longer 
(~ 150 yr). This period of strong mass loss can however not 
be placed closer to the star, since the predicted CO(6-5) line 
becomes far too bright. The estimated mass-loss rate in region 
(D) is ~ 1 x lO^Moyr 1 . 

The estimated kinetic temperature, velocity and mass-loss 
profile yield the predicted CO line profiles displayed in Fig. 8. 
The central peak (with a width of ~ 20 km s~') visible in all 
line profiles originates close to the star (r < 15 R*), where the 
dust is accelerated away due to radiation pressure. Region (A) 
is responsible for the main contribution to the high-excitation 
CO(7-6) and (6-5) lines. With a temperature above ~ 100 K, 
the density in this region results in an optically-thick, non- 
resolved Gaussian line profile (see dashed line in Fig. 9). The 
influence of region (B) is only marginal, and results in a small 
intensity enhancement in all line profiles (see dashed-dotted 
line in Fig. 9). The region with a high density enhancement, de- 
noted as region (C), and with an extent from ~ 600 to ~ 700 R*, 
contributes differently to the high-excitation CO(7-6) and (6- 
5) lines than to the lower-excitation lines. In the former case, 
the temperature (< 90 K) in region (C) is too low to signif- 
icantly populate the higher-excitation energy levels, yielding 
an optically thin, resolved double-peaked line profile, for the 
beamwidths of the JCMT. The situation is markedly different 
for the lower-excitation lines: (1) the low temperature in this 
region causes ~ 5 % of the total CO number density to be in 
the ground state, ~ 13 % to be in the J = 1 state, ~ 21 % in 
the J = 2 state and ~ 23 % in the J = 3 state, yielding optical 
depths larger than 1 . The shape of the line profiles then follows 
the shape of the source function S v , which form can be directly 
inferred from the temperature distribution displayed in the up- 
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per panel of Fig. 7. The increase and subsequent decrease in 
T(r) in region (C) results in a double-peaked line profile, with 
widths in the red and blue peak of some 20 km s _1 . As an illus- 
tration, we have plotted the contribution of region (A)+(B)+(C) 
to the full line profile in Fig. 9. Note that only the (3-2) line 
transition is slightly resolved in region (C), while the larger 
beamwidths for the (2-1) and (1-0) line do not resolve the en- 
velope. It can be seen in Fig. 9 that the intensity in the central 
peak is under-predicted in the CO(2-l) and (3-2) line (dotted 
line in Fig. 9). The low quality of the CO(l-O) line makes it 
difficult to substantiate this for this line. The contribution from 
region (D) fills in this gap, due to the still substantial mass- 
loss rate of M(r) = 1 x 10 _6 Myr _1 in this region and its low 
temperature. 

The small scale ripples visible in the calculated line pro- 
files are not due to numerical noise in the radiative transfer and 
can not be removed by increasing the number of grid-points 
(now taken to 150 gridpoints). The cause for the presence of 
these ripples is the complex temperature structure in region 
(C) which results from the competing effects of heating and 
cooling processes. The associated source function is therefore 
quite complex. In combination with small scale structure in the 
opacity resulting from the number-density profile this explains 
these ripples in the theoretical profiles. However, a less abrupt 
M-profile would cause these ripples to decrease in amplitude. 

The match between observations and theoretical predic- 
tions is not perfect. Main reasons for this are the simplifying 
assumption of a spherical symmetric envelope and a velocity 
structure being continuously increasing with r, as required by 
the comoving frame formalism. In reality the higher mass-loss 
rate in region (C) should drive a dusty wind with a different 
velocity than in the low-mass regions and maybe also with an- 
other dust-to-gas ratio. 

4.5. Comparison with other studies 

A direct comparison with other derived mass-loss rates (see 
Table 3) is difficult, since all of the studies mentioned in Table 3 
assume a density varying as r~ 2 . In case of the analysis of dust 
spectra, one moreover has to assume a dust-to-gas ratio i// to 
estimate the gas mass-loss rate. The main differences with our 
study are that in our case (1) the mass-loss rate is variable, and 
hence the density is not varying as r~ 2 , (2) the dust-to-gas ratio 
is derived self-consistently to obtain the observed expansion 
velocity, and (3) the IR emission traces warm dust close to the 
star, while the analysis of different low-excitation rotational CO 
lines gives us diagnostics on a much broader spatial region. 

Assuming a constant mass-loss rate, the best fit model is 
obtained at M= 3.3 x 10~ 5 M yr~' (see Fig. 10). From the dis- 
crepancy in strength of the higher-excitation rotational CO line 
transitions, it is immediately clear that a model with constant 
mass loss is unable to explain the observed lines. Moreover, this 
model fails completely in fitting the detailed line profiles. This 
is not surprising in view of what is already known about this 
object (see Sect. 4.3). This value for M can be directly com- 
pared to the result of Loup et al. (1993), who derived a value 
of 2.3 X 10~ 5 M yr -1 (scaled to a distance of 1500 pc) from the 



CO(1-0) integrated line intensity. Zuckerman & Dyck (1986) 
used the CO(2-l) integrated line intensity as observed with the 
NRAO 12 m Kitt Peak to derive a much higher mass-loss rate 
of ~ 10~ 4 M o yr _1 . The noise in their observed line profiles 
was so large that the clear composite nature as seen in Fig. 8 is 
not traced. As Zuckerman & Dyck (1986) stated, they scaled 
the analysis of Knapp & Morris (1985) for other stars that 
are located at comparable distances from Earth (~ 1500pc). 
However, Knapp & Morris (1985) derived mass-loss rates from 
CO(1-0) lines of a sample of 50 stars, of which only 2 were 
supergiants. A key point resulting from our study is that the 
analysis of full line profiles with different excitation energies 
provide a much better diagnostics to derive the mass-loss rate 
than peak intensities, also taking into account that the main cal- 
ibration uncertainties do not effect the line shape. 

Values for the mass-loss rate traced by dust character- 
istics are generally larger than our mean mass-loss rate of 
5.1 x 10 _5 M o yr _1 , especially keeping in mind that the dust- 
to-gas ratio used by other authors (see Table 3) is always 
higher than the value we infer from the CO-lines (i// = 0.002). 
However, (near)-infrared dust emission traces regions closer 
to the star than the low-excitation rotational CO lines. The 
mean mass loss rate up to 100R* that we determine is in- 
deed about 0.75(+0.25) x 10" 4 M o yr _1 , substantially higher 
than the mean mass-loss rate averaged over the full envelope 
of 5.1 x lO^Moyr" 1 . 

4.6. Discussion: mass-loss variability 

Smith et al. (2001) also studied mass-loss variability. From 
surface photometry they derived that the star may have ex- 
perienced enhanced mass loss over the past 1000 yr (~ 3 x 
lO" 4 M yr _1 between ~ 0.5" and ~ 7"), with indications of 
a lower mass loss (in the order of 10~ 5 to lO _4 M yr _1 till 
~ 0.5") during the past 100 yr (see third column in Table 2). 
They noted that these estimates can be off by a factor of a few 
due to uncertainties in the calibration of the mid-IR imaging 
data, model assumptions and assumed dust properties. Linking 
our derived mass-loss history with their optical and near in- 
frared images, the enhanced gas density traced in our region 
(A) coincides with several bright nebular condensations located 
at 1" to 2" south and southwest of the star (indicated by 'S' and 
'SW in their Fig. 8). The high mass-loss phase indicated by 
(C) in the bottom panel in Fig. 7 is situated at the same radial 
distance as their 'arc #1', being from 5" to 7" from the star. 
The two other arcs showing up in their Fig. 8 (named as 'arc 
#2' and 'curved nebulous tail') at 2" to 4" from the star, are 
not directly traced in our mass-loss history. However, test com- 
putations show that an extra enhancement of mass loss around 
3" with M < 2 x lO^Moyr 1 with an extent of 1" may be 
present without seriously altering the calculated line profiles 
(see dashed line in Fig. 7, indicated by (E)). Hence, we may 
conclude that our results are consistent with the Hubble and 
near-IR images presented Smith et al. (2001). 

We conclude from our analysis and the study performed by 
Smith et al. (2001), that VY CMa underwent large changes in 
mass-loss rate in the order of a factor 100 — averaged over 
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Figure 9. CO rotational line profiles result- 
ing from the best model (full line) compared 
to theoretical CO line profiles taking into 
account (1) region (A): dashed line, (2) re- 
gions (A)+(B): dashed-dotted line, and (3) 
regions (A)+(B)+(C): dotted line. 
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Figure 10. CO rotational line profiles of V Y 
CMa (Kemper et al. 2003; Nyman et al. 
1992) (grey line) compared with the model 
predictions (black line) based on the stellar 
parameters of VY CMa as given in the sec- 
ond column of Table 2, but with a constant 
M of 3.3 x 10- 5 M o yr' (black line). The 
derived dust-to-gas ratio ip is 0.002. 



all directions — with a time separation of some 800 yr. This 
timeline is of the same order of magnitude as the separations of 
200 - 800 yr observed by Mauron & Huggins (1999) observed 
in the carbon-rich AGB-star IRC +10 216. Also Schoier et al. 
(2005) found that the mass-loss rate producing the faster mov- 
ing wind in some carbon stars with detached shells should be 
almost two orders of magnitude higher than the slower AGB 
wind preceding this violent event. It is important to remark that 
not all dust-driven winds do display this kind of modulations: 
a sample of nine carbon-rich AGB-stars with mass-loss larger 
than 1.5 x 10~ 5 M yr _1 analysed by Schoier et al. (2002) are 
not likely to have experienced any drastic long-term mass-loss 



rate modulations over the past thousands of years. As stated al- 
ready in the introduction, this modulation time scale for AGB 
stars is too large to be due to stellar pulsations, and too short to 
be due to nuclear thermal pulses. It is suggested by Simis et al. 
(2001) that such shells are formed by a hydrodynamical oscil- 
lation in the envelope, while the star is on the AGB. The cause 
of this effect is a subtle mechanism, involving an intricate non- 
linear interplay between gas-grain drift, grain nucleation, radi- 
ation pressure, and envelope hydrodynamics. Another mecha- 
nism proposed by Soker (2000) is a solar-like magnetic activ- 
ity cycle. The periodic change of direction of the stellar mag- 
netic field enables periodic structure formation, with time inter- 
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vals between consecutive ejection events in the order of 200 - 
lOOOyr. In combination with giant convection cells which may 
be present in the outer atmospheric layers of red supergiants 
(Freytag et al. 2002), this may explain the morphology of the 
arcs (Smith et al. 2001). 

5. Summary 

We have developed a code to derive the gas mass-loss history 
from a non-LTE analysis of the full line profiles of CO rota- 
tional transitions with different excitation energies, including 
a self-consistent computation of the velocity structure and gas 
kinetic temperature, taking several heating and cooling mecha- 
nisms in the CSE into account. The code has been benchmarked 
to other existing codes, demonstrating good performance. 

The code is applied to the analysis of the wind structure of 
VY CMa. For the first time, we could pinpoint the mass-loss 
history of this supergiant from the analysis of the CO(1-0), (2- 
1), (3-2), (6-5), and (7-6) line profiles and could demonstrate 
that the mass-loss rate underwent substantial changes during 
the last 1400 yr. We especially note a very strong peak of mass 
loss of ~ 3.2 x 10 _4 M G yr _1 some lOOOyr ago having lasted 
~ lOOyr. This period of high mass loss was preceded by a long 
period of low mass loss, in the order of 1 x 10~ 6 M yr~' taking 
some 800 yr. The current mass loss is estimated to be in the 
order of 5 x 10~ 5 to 1 x lO^Moyr -1 . The mass-loss history 
that we derive is significantly different from what one would 
estimate from a constant mass-loss rate model, and therefore 
paints a much more complete picture on the evolution of the 
mass loss of AGB and red supergiant stars. 
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